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Abstract 

Anisotropic Non-Linear Diffusion is a powerful image processing technique, which allows to 
simultaneously remove the noise and enhance sharp features in two or three dimensional images. 
Anisotropic Diffusion is understood here in the sense of Weickert, meaning that diffusion tensors 
are anisotropic and reflect the local orientation of image features. This is in contrast with the 
non-linear diffusion filter of Perona and Malik, which only involves scalar diffusion coefficients, 
in other words isotropic diffusion tensors. In this paper, we present an anisotropic non-linear 
diffusion technique we implemented in ITK. This technique is based on a recent adaptive scheme 
making the diffusion stable and requiring limited numerical resources. 
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1 Introduction 

A digital image is usually a large two or three dimensional array of pixel values (typically scalars or 
vectors). Image processing methods based on Partial Differential Equations (PDEs) regard images 
as approximations of continuous objects, namely functions from an image domain to a pixel space, 
to which physics-inspired evolution rules can be applied. Among them, Non-Linear Anisotropic 
Diffusion (NLAD) is a variant of the heat equation, generalized in two regards: Non-Linearity and 
Anisotropy . 

Anisotropy in diffusion means that the smoothing induced by the PDE can be favored in some 
directions and prevented in others. This is specified by local eigenvectors and eigenvalues of the 
diffusion tensor field (see §2). Diffusion coefficients are thus location and direction dependent, 
generalizing the approach of Perona and Malik [2] which is only location dependent. Importantly, 
efficient schemes for anisotropic diffusion have been recently made possible by the breakthrough in 
[1]. This has motivated the development of the ITK module presented in this paper. 
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Non-Linearity in diffusion means that diffusion tensors are automatically generated from the 
processed image. We implemented the strategies of Weickert [3] and we give a simple framework 
for designing extensions and variants, see §3. The implemented filters and their parameters are 
described in §4. Figures 2 and 3 illustrate their effect on 2D images; Figures 4 and 5 on 3D images; 
Figures 6 and 7 on color and vector images. 

A possible application of NLAD is to enhance a fingerprint image by smoothing tangentially 
to the lines. Evidence is also plentiful for NLAD relevance in many other image processing ap¬ 
plications, but its use has been limited by technical aspects so far. We intend to alleviate such 
limitations with the present contribution. 

Notations: Let d G {2, 3} denote the image dimension, let fl C be the image domain, and let 
V be the pixel space (e.g. [0,1] for grayscale, [0, l] 3 for color, R d for vectors). Throughout the paper, 
we informally consider an idealized cartoon image model, involving a set V C D of image contours 
of dimension d — 1. The processed image u : D — >> V is smooth on O \ T, but has discontinuities 
across T, and is overall corrupted by e.g. by additive white noise. A key feature of NLAD is its 
ability to detect the set T and smoothen tangentially to it. Finally, let S'j" denote the collection of 
symmetric positive definite d x d matrices, and let Id be the identity matrix. To each D G S'j" we 
associate the norm \\e\\]j := yj (e, De), eGR d , where (*, •} denotes the standard scalar product on 
R d . 


2 Linear Anisotropic diffusion 

Linear Anisotropic Diffusion 1 (LAD), in divergence form, is an elliptic PDE which reads 

dtu = div(DV'u), (1) 

where D : D S j" is a given field of symmetric positive definite diffusion tensors. Eigenvectors 
of these tensors define preferential diffusion directions, and the eigenvalues their corresponding 
coefficients. Evolution rule (1) is complemented with an initial condition u( 0, •) = uq at time 
t — 0. If uq has pixels of vector type, then their components are treated independently. We use 
Neumann 2 conditions on the domain boundary dD, as is common in image processing. LAD is 
formally a continuous gradient descent for the elliptic energy 

£{u) := f ||Vu(x)||^ (a!) dx. (2) 

J Cl 

Qualitative effects of LAD strongly depend on the chosen field D of diffusion tensors. Choosing D = 
Id identically on D yields the standard heat equation, which qualitative properties are well known 
in image analysis: any noise present in the image u is quickly eliminated, but in the meanwhile all 
image sharp features are blurred. 

This undesirable side effect can be limited with a proper choice of diffusion tensors D : D — >► 
Sj . Indeed LAD smoothes primarily the image features which contribute strongly to the energy 
(2). In the spirit of Perona and Malik, one can introduce an isotropic but variable conductivity 

D (x) = c(x) Id, with c(x) <C 1 for all x close to the image contours T, see introduction. Smoothing 

is prevented in the neighborhood of T, which preserves the contours sharpness, but also traps some 

1 We use here the terminology of Weickert [3]. Perona-Malik diffusion, which uses an adaptive scalar tensor field 
similar to (4), is in contrast a Non-Linear Isotropic Diffusion equation. 

2 Neuman boundary conditions must take into account the geometry defined by the diffusion tensor field. They 

take the form (Vw(a:),D(x)n(x)) = 0, where n(x) denotes the unit outward normal at x G <9P. 
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noise along them, see Figure 2 (IV). A more elaborate approach is to construct anisotropic diffu¬ 
sion tensors D(x), which favor diffusion tangentially to the contours curves T, but simultaneously 
prevent diffusion transversally to these curves and between different image regions. All noise is 
eliminated, yet image discontinuities are preserved, see Figure 2 (II). 

Numerical schemes for LAD are in general non-trivial due to interaction between the anisotropic 
geometry of the diffusion tensors, and the cartesian structure of the pixel grid. The authors recently 
developed [1] a numerical scheme which handles this interaction using special tools from discrete 
geometry, named Lattice Basis Reduction (LBR). It provides strong mathematical guarantees (con¬ 
sistency, stability, maximum-principle) for a limited numerical cost. 

3 Non-Linear Anisotropic Diffusion. 

Linear Anisotropic Diffusion, discussed in §2, requires two main inputs: an image uo serving as 
an initial condition, and a field of diffusion tensors. In order to reduce user input, the diffusion 
tensors D can be defined in terms of the filtered image u. The resulting PDE is called non-linear 
anisotropic diffusion 

0 t u = div(D w Vi/), (3) 

complemented, again, with Neumann boundary conditions. Perona and Malik [2] suggested to use 
the following non-linear isotropic (i.e. proportional to the identity matrix) tensors 

t) := c u (x, t) Id, with c u (x,t) :=— . =, (4) 

V ' V 7 y/l + \\Vu(x,t)\\ 2 /\*' ^ ' 

where A > 0 is a user specified constant. Diffusion is prevented where the conductivity c u (x,t) 
is small, in other words where \\Vu(x,t)\\ is large, such as along the image contours T. Perona- 
Malik diffusion is already available 3 in ITK. It has been the subject of considerable academic and 
industrial interest revealing that, in spite of its numerous qualities, it is mathematically ill posed, 
unstable, often leads to unsightly “staircasing” visual artifacts, and is not adequate for oscillating 
patterns as in Figure 3. 

We describe in the following Coherence Enhancing Diffusion (CED) and Edge Enhancing Dif¬ 
fusion (EED), which are based on more complex tensor constructions introduced by Weickert [3]. 
Our first ingredient is the Gaussian convolution kernel: given a standard deviation a > 0 

K *( x ) : =jd K i0)’ where K ^ x ) : = (^y exp (^^) • ( 5 ) 

The structure tensor S u : ft S'j", defined below, is a robust estimator of the gradient direction in 
an image u , even if this image has oscillating textures. It depends on two small positive parameters: 
the noise scale cr, and the feature scale p. We denote by * the convolution operator, and by 
v ® v = vv T the self outer product, which yields a semi-definite symmetric matrix. 

S u := K p * (V Uq- ® V?v), where u a := K a * u. (6) 

If u is an image with vector pixels, then S u is the sum of the structure tensors associated to the 
components of u. Assume that u is a scalar image, fix a time t and a point x G D, and denote 
S := S u (x, £), v := Vu(x, t). Let also Ai < • • • < denote the eigenvalues of 5, sorted by increasing 
magnitude, and ei, • • • , the corresponding unit eigenvectors. If u is sufficiently smooth, then the 

3 Under the name (misleading with our conventions) Gradient AnisotropicDiffusionlmageFilter . 
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largest eigenvalue approximates the gradient squared norm: A^ ~ ||p|| 2 , while the corresponding 
eigenvector approximates the gradient direction: ~ 

Weickert’s diffusion tensors D := ~D u (x,t), are defined in terms of this eigen-analysis of the 
structure tensor S := S u (x,t): 

D = ® e ii where S = ^ A*e* ® e*. (7) 

1 <i<d 1 <i<d 

Smoothing is promoted in the direction if ^ is large, and prevented if fii is small, for any 
1 < i < d. Weickert’s classical constructions are presented in (8) and (11). One should not shy 
away of designing more complex and application dependent variants; for instance one may want to 
enhance filaments and tubular structures in 3D data. Three very simple variants (9), (10) and (12) 
are presented for illustration. All depend on three parameters A, m, a. The main one, A > 0, is an 
edge detection threshold. The exponent m is typically 2 or 4. The small parameter cp typically 
1/100, determines the condition number of the diffusion tensors. 

• Edge Enhancing Diffusion (EED) aims to avoid significant diffusion across the set T of image 
contours, but to allow it anywhere else. Note that for 3D images discontinuity planes will 
be enhanced, rather than edges. The first diffusion tensor eigenvalue is /ii = 1, because the 
eigenvector c\ is orthogonal to the image (approximate) gradient direction e^, hence never 
transverse to T. Other eigenvalues satisfy fii « a <C 1 if A^ — Ai > A, and ^ « 1 otherwise. 
The condition A^ — Ai > A indeed suggests that the eigenvector points through an image 
contour. Precisely 4 : (note that \±\ — 1) 

Mi := 1 - (1 - a) exp (- ^ _ A x ) ) ' ( 8 ) 

The choice of Weickert, to set fi\ — 1, may lead to undesired effects: one always performs 
diffusion in at least one direction. An undesirable side effect is that the image is blurred close 
to the angles of its contour set T. We believe that such salient features should be preserved, 
hence we introduce a Conservative variant of EED (cEED) for which {i\ can be small, when 
appropriate, so as to prevent diffusion around the angles of T, see Figure 8 for a comparison. 
Precisely 

exp (- (£) 

If all eigenvalues are set equal fii = • • • = /i^, then the diffusion tensors are isotropic , in other 
words scalar multiples of the identity. The following isotropic variant of EED is close in spirit 
to the Perona-Malik model: diffusion is prevented in the neighborhood of the image contours 
r, regardless of direction. This construction is implemented purely for comparison with the 
anisotropic ones, and does not take advantage of the innovative numerical scheme developed 
by the authors 

IM := 1 - (1 - a ) exp (A) ) . (10) 

• Coherence Enhancing Diffusion (CED) prevents diffusion except along local image structures 
which have a coherent direction. The diffusion tensor eigenvalues satisfy « a 1, unless 

4 Actually, Weickert uses pi = 1 together with (9) for i > 1. This results in discontinuous diffusion tensors, which 
is not advisable from a mathematical standpoint, hence the formula (8). 


Hi := 1 - (1 - a) 
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// Usage with uchar pixel type, cast to double, 
typedef Imagecunsigned char,3> ImageType; 
typedef Image<double,3> FPImageType; 

typedef CastImageFilter<ImageType, FPIinageType> CastFilterType; 
typedef CoherenceEnhancingDiffusionFilter<FPImageType> CEDFilterType; 

CastFilterType: : Pointer castFilter = CastFilterType: :New(); 
CEDFilterType: : Pointer cedFilter = CEDFilterType: :New(); 
cedFilter->Set!nput(castFilter->GetOutput()); 


// Usage with RGB pixel type, cast to Vector 
typedef Vecto refloat, 3> VectorType; 
typedef RGBPixel<unsigned char> RGBType; 
struct CastFunctorType { 

VectorType operator () (const RGBType & x){ 

VectorType result; for(int i=0; i<3; ++i) result [i] =x[i]; return result;} 

}; 

typedef ImageeRGBType, 2> ImageType; 
typedef ImagecVectorType, 2> FPImageType; 

typedef UnaryFunctorImageFilter<ImageType, FPImageType, CastFujnctgrType> CastFilterType; 
typedef CoherenceEnhancingDiffusionFilter<FPImageType, ftoat> CEDFilterType; 


Figure 1: The provided filters expect image pixels of scalar or vector pixel type. A preliminary 
cast is required for integral pixel types (left), or non-Vector pixel types such as RGB (right). The 
provided executable (or “.cxx” file) performs these casts automatically. 


if Ad — A^ > A in which case f±i « 1. The condition A^ — A^ > A indeed suggests that points 
through an image feature, and that points tangentially to it. Precisely: (note that ^ = ck) 

m:=a + (l- a) exp (- ^ ^ ^ ) ). (11) 

The above formula often leads to false positives: at a position with large gradients, but 
without a clear preferred direction, one may very well have A^ — A^ A (for instance if 
A d = 100A and A^ = 95A). A more reliable coherence detector is A m — A^ >> A + A^, which 
leads to a Conservative variant of CED (cCED), see Figure 8 for a comparison. Precisely: 
(note that fid = a) 

/Xi :=a + (l-a)exp(-^T^T) ). (12) 

We emphasize that the distinction between EED and its conservative variant cEED (resp. CED 
and cCED) is rather subtle, and mostly located around image contour corners, as evidenced on 
Figure 8. In other illustrations, we only show the conservative variant, which is slightly better at 
preserving detail. 

4 Implemented filters and their parameters 

Our contribution to ITK consists of the following four image filters, which implement the mathemat¬ 
ical notions presented in §2 and §3. The figures are produced with the last filter, which implements 
CED, EED and their variants described §3. The first three filters are its building blocks, but they 
may be of independent interest for other applications. All filters are multithreaded. 

For each filter, the processed image pixels can be of scalar or vector type. In the latter case, the 
underlying floating point type needs to be specified via the second template parameter of the filter, 
see Figure 1 (bottom right). Pixels of integral type (resp. RGB pixels) must be cast to floating 
point (resp. Vector ) types, see Figure 1. The image dimension must be 2 or 3. 

Linear Anisotropic Diffusion (LAD). The filter LinearAnisotropicDiffusionLBRImage- 
Filter , requires two inputs: a processed image and a tensor image. Note that non-linear 
diffusion is achieved through successive linear diffusions, over multiple small time intervals, 
with regularly updated diffusion tensors. Parameters: 

• MaxDiffusionTime specifies the target physical time for the LAD evolution PDE (1). 
The filter has early abort options, as discussed in the next point, hence one should check 
the EffectiveDiffusionTime at termination. 
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• MaxNumberOfTimeSteps , RatioToMaxStableTimeStep. Explicit numerical schemes for 
diffusion are subject to a Courant-Friedrichs-Levy (CFL) condition, which limits the 
largest stable time step. This time step, which depends on the input diffusion tensors, is 
automatically computed by the filter, and the requested time interval [0, MaxDiffusionTime\ 
is split accordingly. Early abort occurs if this splitting exceeds the specified MaxNum¬ 
berOfTimeSteps. 

Structure tensor. Filter StructureTensorlmageFilter . Parameters: 

• NoiseScale cr, and FeatureScale p , see (6). Suggested ranges: a G [0.5,3], A G [2,10], 
assuming a unit pixel spacing. The lower bounds of these intervals are recommended, 
unless noise is extremely strong. 

• RescaleForUnitMaximumTrace. If on, the structure tensors are rescaled: S' u = aS u , 
where a > 0 is the largest constant such that Trace(S^(x)) < 1 for all x G D. (The 
trace is the sum of the eigenvalues Ai + • • • + Asee (7).) This option is meant to ease 
the choice of the edge detection threshold A in CED, EED. One may want to check the 
variable PostRescaling = a at termination. 

Non-Linear Anisotropic Diffusion (NLAD). Filter AnisotropicDiffusionLBRImageFilter 

. This class (I) computes structure tensors by invoking the previous filter, (II) performs their 
eigen-analysis, (III) changes them into diffusion tensors via (7), (IV) runs linear diffusion 
with the constructed tensors by invoking the first filter. After a limited number of time-steps 
of linear diffusion, the steps (I-II-III-IV) are repeated so as to update the diffusion tensors, 
until exhaustion of the prescribed diffusion time. Parameters: 

• DiffusionTime for which the evolution rule (3) is applied. 

• Adimensionize. If on, the filter ignores the image pixel spacing information, and sets 
on the RescaleForMaximumUnitTrace option for structure tensor generation. This is 
intended to ease the choice of DiffusionTime , and of the edge detection threshold A in 
CED, EED, and variants. 

• MaxTimeStepsBetweenTensorUpdates is self descriptive. NoiseScale and FeatureScale 
are passed for structure tensor generation. 

• Eigenvalues Transform is a virtual method used to construct the diffusion tensor eigen¬ 
values from those of the structure tensors (Awhich are sorted increasingly 

for convenience. The method must be redefined in a subclass, as in the next filter, else 
it triggers an exception. 

Coherence-Enhancing diffusion and Edge-Enhancing diffusion. The two PDEs, and their 
variants, are implemented in the same filter CoherenceEnhancingDiffusionFilter , which 
subclasses the Non-Linear Anisotropic Diffusion filter. Parameters: 

• Enhancement allows to switch between EED, cEED, CED, cCED and Isotropic (10) 
tensor constructions, by redefining the superclass virtual method EigenValuesTransform. 
The relevant choice depends on the type of image structures that one wants to enhance. 

• Lambda— A, Exponent= m, Alpha= a are the parameters involved in tensor design (8) 
-(12). 

• Suggested parameter ranges. Adimensionize— True. DiffusionTime G [0.5, 5], although 
larger values can be relevant for very strong noise or artistic effects. Edge detection 
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Figure 2: I: Source image. II: cEED is the best pick. Ill: cCED only affects the directionally 
coherent parts, the image region contours, which is clearly inadequate. IV: Isotropic (10, Perona- 
Malik like) diffusion leaves some noise trapped along the image contours. Parameters: T = 20, 
A = 0.05, a — 3, others default. 



Figure 3: I: Source image. II: cEED is not advised, since it blurs the junctions of the fingerprint 
lines. Ill: cCED is the best pick, since it enhances the fingerprint lines, and does not blur the 
more complex regions. IV: Isotropic (10, Perona-Malik like) diffusion either fails to remove noise, 
or blurs the fingerprint lines, depending on the image region. Parameters: T — 20, A = 0.02, others 
default. 


threshold A G [10 -3 , 5 x 10 -2 ], small for complex images with detail, large for sim¬ 
ple “cartoon” like images. Finally m G [2,4], a — 0.01, though these parameters are 
secondary and have little impact. 
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Figure 4: Left: Volume plot and slice of an FMRI of the human skull, with artificially added 
Gaussian noise of variance 0.01 (data range: [0,1]). In the volume plot, values below 0.08 are 
shown transparent. Right: Effect of cEED. Parameters: T = 5, A = 0.003. 



Figure 5: Synthetic function cos(||x|| 3 ), x E [0, l] 3 , corrupted with Gaussian noise of variance 0.5. 
Left: level lines, and slice {x = 0.1}. Right: effect of CED. By construction, CED barely affects 
the neighborhood of the origin (0,0,0), where no specific coherent direction can be determined. 
Parameters: T = 10, A = 0.02, a — 4, p — 10. 



Figure 6: I: Detail of the Lena image. II: cEED removes noise and preserves most image detail. Ill: 
cCED enhances the flow structure of the hat plumes, without much affecting the rest of the image. 
In particular it does not remove noise. IV: Isotropic diffusion (10, Perona Malik like) removes noise 
in the interior of the image regions, but leaves some noise close to image contours such as the border 
of Lena’s cheek. Parameters: T — 2, A = 0.003, m — 4, others default. 
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Figure 7: Left: Directions and norms (black=0, white=l) of the vector field v(x) := sign(||x|| — 
1) ar L /||:c|| on [—1.3,1.3] 2 , degraded by gaussian noise of variance 2, on a 50 x 50 grid. Norms of 
these vectors (black=0, white=l). Right: Effect of cEED. Note that the streamlines are better 
reconstructed, and that the norms vanish along the (approximate) circle ||x|| = 1 where v(x) 
changes sign so that there is a cancellation effect (and likewise close to the vector field singularity 
at the center). Parameters: T — 10, A = 0.05, others default. 



Figure 8: Left: (I) source image, (II) effect of cEED, (III) effect of EED. Note that the triangle 
corners are better preserved with cEED. Parameters: T — 5, A = 0.03. Right: (I) source image, 
(II) effect of cCED, (III) effect of CED. Note that cCED better preserves contrast along where the 
two fronts meet. Parameters: T — 20, A = 0.05. 
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